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A self-consistent global fitting method based on the Markov Chain Monte Carlo 
technique to study the dark matter (DM) property associated with the cosmic ray 
electron/positron excesses was developed in our previous work. In this work we 
further improve the previous study to include the hadronic branching ratio of DM 
annihilation/decay. The PAMELA p/p data are employed to constrain the hadronic 
branching ratio. We find that the 95% (2a) upper limits of the quark branching 
ratio allowed by the PAMELA p/p data is ~ 0.032 for DM annihilation and ~ 
0.044 for DM decay respectively. This result shows that the DM coupling to pure 
leptons is indeed favored by the current data. Based on the global fitting results, we 
further study the neutrino emission from DM in the Galactic center. Our predicted 
neutrino flux is some smaller than previous works since the constraint from 7-rays 
is involved. However, it is still capable to be detected by the forth-coming neutrino 
detector such as IceCube. The improved points of the present study compared 
with previous works include: 1) the DM parameters, both the particle physical 
ones and astrophysical ones, are derived in a global fitting way, 2) constraints from 
various species of data sets, including 7-rays and antiprotons are included, and 3) 
the expectation of neutrino emission is fully self-consistent. 
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I. INTRODUCTION 



In a previous work Paper I) we have developed a Markov Chain Monte Carlo 
(MCMC) code to fit the parameters of dark matter (DM) models, which are proposed to 
explain the recent reported abnormal excesses of cosmic ray (CR) positrons and electrons by 
PAMELA Q, ATIC Q, HESS 5| and Fermi-LAT 6]. One assumption adopted in that 
work is that the DM particles only couple with leptons. To extent the discussion including 
hadronic channels will be natural and necessary. Actually the non-excess of PAMELA p/p 
data 131 have been studied in many works to set constraints on the hadronic couplings of DM 

nnn 

particles (e.g., [8l. 191. llOj). However, these studies can only give some illustration constraints 



on the model parameters instead of a full scan of the possible parameter space. Our MCMC 
fitting scheme makes it possible to fully scan the high-dimensional parameter space and 
derive the model-independent constraints from the data directly. 

Shortly after the proposal of using DM annihilation to account for the data, the accom- 
panied 7-ray and radio emission are investigated as constraint and/or future probe of the 

nnn 

current models (e.g., LLU [121 ll3||). Besides the photon emission, neutrino can also be re- 

nnn 

garded as an another probe to test the DM interpretation of the data [14], Ufa, H6[ • In Ref . 



141 ] it is shown that Antares and IceCube are promising to detect the neutrino signals of 



DM annihilation in the Galactic center (GC) or the massive subhalo for models to explain 
PAMELA and ATIC data, assuming an NFW profile 17| of DM spatial distribution. It is 
also claimed that for annihilation DM scenario the IceCube/DeepCore detector can explore 
much of the parameter space to explain the PAMELA and Fermi-LAT e ± data, still for an 
NFW profile 15|. For decaying DM scenario IceCube/DeepCore also has the potential to 
exclude some of the parameter space, depending on the final states and DM profile [la ]. In 
all of these works, the expectation of neutrino signals will strongly depend on the annihila- 
tion/decay final states and DM profile, and the constraints from various kinds of data are 
not taken into account simultaneously. 

Based on the above points, in this work we improve our MCMC code to include the 
hadronic branching ratio, and then we re-examine the neutrino signals according to the 
parameter sets derived in MCMC calculations. We include the positron fraction data from 
PAMELA the electron + positron data from Fermi-LAT and HESS {4, 5, 6], the diffuse 
7-ray data of the GC ridge from HESS [lsl . and the antiproton-proton ratio data from 
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PAMELA Q in this MCMC study. 

The paper is organized as follows. In Sec. II we introduce the production and propagation 
of the CR electrons/positrons and antiprotons. In Sec. Ill we give the MCMC fitting results 
of DM model parameters, for both annihilation and decaying scenarios. The neutrino signals 
from GC are discussed in Sec. IV. Finally, Sec. V is the summary. 



II. PRODUCTION AND PROPAGATION OF e ± AND p 
A. Solution of propagation equation 

In the Galaxy the transport of charged particles is affected by several processes. The 
scattering off random magnetic fields will lead to spatial and energy diffusions. The stellar 
wind may also blow away the CRs from the Galactic plane. In addition, interactions of CR 
particles with the interstellar radiation field (ISRF) and/or the interstellar medium (ISM) 
can result in continuous and catastrophic energy losses. Since the detailed processes affect 
the propagation are species-dependent, the treatments for positrons and antiprotons are 
separated. For the transport processes we take a spatial independent diffusion coefficient 
D(E) = (3DqR, 5 (where 1Z = pc/Ze is the rigidity) and a constant wind V c directed outwards 
along z. CRs are confined within a cylinder halo L, i.e. the differential density is bound by 
n(z = ±L, -Rmax) = with -R max the scale of the visible Galaxy. The free parameters of the 
model are the halo size L of the Galaxy, the normalization of the diffusion coefficient D 
and its slope 5, and the constant galactic wind V c . 

The propagation equation of CRs can be generally written as 

f) JV f) ( d W \ 

- DAN + V C — + 2hT tot 5(z)N + — —N = q(x, E), (1) 



dz \ / dE y df 

where r tot = ^2 i=H He n-i Oi v is the destruction rate of CRs through interaction with ISM in 
the thin gas disk with half height h ~ 0.1 kpc, dE/dt is the energy loss rate, and g(x, E) is the 
source function. The solutions for electrons/positrons and antiprotons are presented in the 
Appendix. More details about the propag ation processes and the solutions of propagation 



equation can be found in Refs 



20 



21 



22 



23]. 



For the propagation parameters, we use the medium (referred as "MED") set o 
eters which is derived through fitting the observational B/C data given in Ref. 



param- 



m, i.e., 



i 



D = 0.0112 kpc 2 Myr" 1 , 5 = 0.70, V c =12 km s" 1 and the height of the diffusive halo L = 4 
kpc. 



B. Background 

For all the CRs we consider here, there are backgrounds originated from the traditional 
astrophysical sources and/or interactions in the Milky Way (MW) or the Earth atmosphere. 



The CR proton and 
rameterizations in Ref. 



ielium spectra are well measured at the Earth. We adopt the pa- 
s| as the interstellar proton and Helium spectra. Then we calculate 
the positrons (together with a secondary electron component with almost the same flux) 
and antiprotons produced through interactions between CR protons, Helium and the ISM. 
The interaction is restricted in a thin disk with half height ~ 0.1 kpc and the average ISM 



density is adopted as 



best-fit to the B/C data 



. cm 3 . These parameter choices were shown to be able to give 



251 ] . The shapes of secondary positrons and antiprotons are fixed 
to the calculated results, and we further employ two normalization parameters c e + and Cp to 
describe the uncertainties about the inelastic hadronic cross section and propagation effect. 



Cp is found to be within 1 ± 0.25 [8j. For c e + we restrict it in a larger range of 1 ± 0.5 since 



the positrons are more sensitively dependent with the propagation 26|. 

The background of primary electrons is different from that of positrons and antiprotons. 
Since we do not have exact knowledge of the primary electron injection spectrum from the 
acceleration source, we adopt a 2-parameter power- law function q e - = a e -E e ^ e ~ to describe 
the injection source of primary electrons. This is similar to the case of GC 7-rays, which are 
parameterized by (f)^ kg = a 1 E 1 b "'. a e -, b e -, a 7 and b 1 are fitted in the MCMC procedure. 

There should be two kinds of backgrounds of neutrinos: the atmospheric background 
and the astrophysical one. For energies smaller ~ 100 TeV the atmospheric background is 



dominant |27|]. Thus we only consider the comparison with atmospheric background in this 
work. The result of atmospheric neutrino background is adopted as the direction-dependent 
calculation 



28 



based on the muon data. 
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C. DM contribution 

The source function (emissivity) of DM annihilation or decay to standard model particles 
can be written as 

(CT '" >dAr Vw, (2) 



' 2mi ah 

i X- 



for DM annihilation, or 

" P(r), (3) 



for DM decay, where m x is the mass of DM particle, (crv ) or r is the annihilation cross 
section or decay age of DM respectively, B{ is the branching ratio to final state channel i, 
is the yield spectrum of j species for one annihilation or decay for channel i, and p(r) 
is DM spatial density in the MW halo. In this work we consider three channels to lepton 
pairs e + e~, r + r~ as well as a quark channel qq. Since the antiproton production 

spectra from decay of various quark flavors do not differ much from each other, we do not 
distinguish quark flavors but use an average result from all the flavors. The spectra ^|. 
of positrons, antiprotons, 7-photons and neutrinos from decay of the final state particles of 



various channels are calculated using PYTHIA package 29]. 



Similar as in Paper I, we take the density profile of the MW halo as the form 

p(r) = ^ U) 

where 7 represents the central cusp slope of the density profile, r s and p s are scale radius 
and density respectively. For the MW DM halo, we adopt the total mass to be Mmw ~ 10 12 
30 ] and the concentration parameter to be cmw ~ 13.5 3l|. Then we have r s = 



'3 



?"mw/cmw(2 — 7) where tmw ~ 260 kpc is the virial radius of MW halo. Then p s can be 
derived by requiring Mmw = / pd.V. The local density p in this process is checked to be 
within 0.27 to 0.25 GeV cm -3 for 7 varying from to 1.5. 

Finally there are propagation effects of these particles before being detected by CR de- 
tectors. For charged particles such as positrons and antiprotons, we can replace the source 
term in Eq. ([1]) with Eqs. ([2]) ([3]) and solve the propagation equations to get the propagated 
fluxes. For neutral particles like 7-ray photons and neutrinos we just need to integrate the 
contribution along the line-of-sight of given direction. More details of the treatment of 7- 
rays can refer to Paper I. Note that for neutrinos there will be oscillations between different 
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flavors. We actually count all flavors of neutrinos from the output of PYTHIA, and multiply 
a factor 1/3 to give the result of muon neutrinos. 

III. MCMC STUDY OF DM PARAMETERS 

The full parameter space of this MCMC study is 

P = (m x , (av) or r, B e , B T , B q , 7, a 7 , 6 7 , a e -, b e -, c e +, c p ). (5) 

For the branching ratios we have a normalization condition B e + B^ + B T + B q = 1. As a 
consequence we have 12 free parameters in total. 

The probability distribution of the fitting parameters is shown in Fig. [TJ The basic fitting 
results are similar with our previous study in Paper I. What's new is the constraint on quark 
branching ratio from PAMELA p/p data. We find that the 2a upper limits of B q is 0.032 
for DM annihilation and 0.044 for DM decay respectively. These results show that the DM 
to explain the recent CR data indeed needs to dominantly couple with leptons instead of 
hadrons. 




FIG. 1: Probability distributions of the model parameters in annihilation (left) and decaying (right) 
DM scenarios respectively. 
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IV. NEUTRINO EMISSION 



In this section we discuss the neutrino emission of the DM models. For each parameter 
set in the MCMC samples obtained in the previous calculation, we compute the neutrino 
flux as a function of energy. The probability distribution of neutrino flux then can be derived 
for each energy bin. 

The predicted fluxes of + from the GC direction are shown in Fig. [2j Two sky 
regions, 1° and 60° (half angle of the cone) around the GC are calculated to show the effects 
of angular resolution. For comparison we also show the atmospheric background adopted 
from Ref . 28( . The average results of atmospheric neutrino fluxes are almost the same for 
1° and 60° cones. 





FIG. 2: Predicted fluxes and la uncertainties of muon and anti-muon neutrinos from the GC for 
annihilation (left) and decaying (right) DM scenarios respectively. The blue solid line shows the 
atmospheric background. 



It is shown that generally the neutrino fluxes are dominated by the atmospheric back- 
ground. For the decaying DM scenario the results are almost always lower than the atmo- 
spheric background, even for the angular resolution as good as 1°, which corresponds to the 
best performance of IceCube 32|. For the annihilation DM case, the neutrino flux from DM 
will have a chance to exceed the atmospheric background for energies ~TeV, given a very 
good angular resolution. Otherwise the signal will be still dominated by the background. 

However, to better see the detectability of such DM-induced neutrinos, we need to com- 



s 



pare the muon events induced by the neutrinos on a detector. Following the method given 



in Ref. [33j we calculate the expect number of muon events on IceCube. Both the con- 
tained and through-going events are included. The size of IceCube detector is adopted as 
~ 1 km 3 , and the data accumulative time is taken as 3 years. The expected muon events 
from both the atmospheric and DM-induced muon neutrinos are compiled in Table [H In 
the calculation we use the mean values of neutrino fluxes shown in Fig. [2j The integral 
energy for muons is adopted from 0.5 to 2 TeV, which ensures us to include DM-induced 
signal and exclude atmospheric background as effective as possible. It is shown that for a 
circle sky region with radius 1° the statistical significance of DM signal is less than 5a for 
IceCube running for 3 years. The case of DM annihilation is some better than DM decay. 
If we enlarge the sky region to about 60° (ir sr), the number of events will be much larger, 
both for atmospheric background and DM signal. The significance of DM signals can reach 
~ 25 — 35a. However, we should keep in mind that in this case it would be still difficult to 
pick out signal events from the background events. The information of spatial distribution 
and energy spectroscopy is necessary. 

TABLE I: Number of muon events between 0.5 and 2 TeV induced by f„ + from DM and 
atmospheric background, for 3-year data taking of IceCube. The number in parentheses is the 



detection significance defined as iVDM/V^i 



cone angle atm DM-ann(<r) DM-decay(cr) 
1° 23.4 20.4(4.2) 6.1(1.3) 

60° 7.68 x 10 4 9.71 x 10 3 (35.0) 7.07 x 10 3 (25.5) 

Compared with the results given in Ref. [jjj], our predicted neutrino flux for DM anni- 
hilation is several times smaller 1 . This is because in our work the 7-ray constraint is taken 
into account and a smaller value of the slope of DM central cusp 7 is required. Whereas in 
Ref. [14I NFW profile is adopted. For decaying DM scenario the results do not differ much 
from each other. Our prediction of neutrino fluxes should be more self-consistent and more 
reliable. 



1 Note this is a rough comparison because the model parameters in Ref. [3] are different from our global 
fitting ones. 
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V. SUMMARY 

Based on the MCMC code developed in Paper I, we further study the properties of 
DM models which may connect with the recent CR lepton excesses, after incorporating the 
hadronic branching ratio and the constraint from PAMELA p/p data. Our global fitting 
results show that DM with annihilation or decay channels dominantly to the combination of 
r + r~ and fi + fi~ can well describe the electron and/or positron data measured by PAMELA, 
Fermi-LAT and HESS. The p/p data from PAMELA limit the branching ratio to quark 



pairs to i 
studies 



)C 



ess than 5%. This conclusion is qualitatively consistent with the previous 



js,^ loj |. but should be regarded as the first quantitative result based on the global 
fitting method. 

We then calculate the neutrino fluxes from the DM annihilation or decay according to 
the MCMC fitting parameters. Due to the 7-ray constraint on the DM profile, Our expected 
neutrino flux for DM annihilation scenario is smaller than that in the previous works where 



the NFW profile is priorly adopted. However, as we have shown in Ref. [34j and Paper 
I, the NFW profile will be strongly constrained by the diffuse 7-rays from the GC ridge 
observed by HESS [JiJ for DM annihilation scenario. The derived 2a upper limit of the 
slope of DM central profile is ~ 0.5 in the MCMC study, which is much smaller than the 
NFW profile with 7 = 1. Although smaller, the neutrino flux still could be detected by 
the forth-coming detector such as IceCube. However, a careful analysis about the angular 
distribution and spectral distribution is necessary to separate the DM-induced signal from 
the high atmospheric background. For DM decay scenario the constraint from 7-rays is 
weaker, and the differences of the neutrino flux between this work and other works are 
smaller. 

APPENDIX: SOLUTIONS OF THE PROPAGATION EQUATIONS OF 
POSITRONS AND ANTIPROTONS 

1. Positrons 

The dominant process in the propagation of positrons is energy loss due to synchrotron 
and inverse Compton scattering for energies higher than ~GeV. In this paper we will neglect 
the convection and reacceleration of positrons, which is shown to be of little effect for E > 10 



10 



GeV (also the interested energy range here) 26|. Then the propagation equation is 

- DAW + s§(f") = «fr^ < A1 > 

in which the second term in the left hand side represents the energy losses. The energy 
loss rate of positrons due to synchrotron and inverse Compton scattering in the MW can 
be adopted as dE/dt = -e 2 /r E , with e = E/l GeV and r E « 10 16 s |35J]. We directly 
write down the propagator for a point source located at (r, z) from the solar location with 



monochromatic injection energy Eg 



22 



23] 



ge & + (r,z,E^E s ) = ^-xg Q (r,z,i) 



in which we define a pseudo time f as 



j-i s-i 



T = T E - 



1-5 



(A.2) 



(A.3) 



& (r, z, f) is the Green's function for the re-arranged diffusion equation with respect to the 
pseudo time f 

' xg 1D (z,f). (A.4) 



A , - 0(t) 

G®\r,z, t) = ; exp 

The effect of boundaries along z = ±L appears in (? 1D only. Following Ref. 
distinct regimes to approach Q m : 



221 ] we use two 



for C, = L /AD t 1 (the extension of electron sphere A = \JADqt is small) 



g m ( z ,f)= J2 (-!) n 



VAttD^t 



exp 



4A)T 



(A.5) 



where z n = 2Ln + (—l) n z; 



otherwise 
G m (z,f) 

where 



- ^ [ex V {-D Q klf)<j> n {Q)<j> n {z) + exp(- J D A;; 2 r)^(0)^(z)] , (A.6) 



n=l 



i> n (z) = sm[k n (L - \z\)}; k n = (n - l/2)n/L, 
VJz) = sin[k' n (L - z)}; k' n = nn/L. 



(A.7) 
(A.8) 



For any source function q(r, z, 8; Eg) the local observed flux of positrons can be written as 



$0 = — x 2 / dz 



Ait 



rdr / dEsGgir, z,E^- E s ) / d9q(r, z, 9; E s ). (A.9) 



2?r 
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2. Antiprotons 

It has been shown that for the propagation of antiprotons neglecting the continuous 
energy losses and reacceleration can provide a good enough approach, especially for energies 



higher than several GeV 20|. We will also adopt this approximation here. Therefore the 
relevant processes include the diffusion, convection and the catastrophic losses — inelastic 
scattering and annihilation in interactions. The propagation equation is 

- DAN + V c ^ + 2hT tot 5(z)N = q(x, E), (A.10) 

where r tot = ^2 i=H He ni o f v is the destruction rate of antiprotons in the thin gas disk with 



half height h ~ 0.1 kpc 20J, g(x, E) is the source function. The propagator for a point 



source located at Xs, expressed in cylindrical coordinates (r, z) (symmetric in 8) is 20] 



S P e(r,z,E) = ^Jff x f^^o (rVW+M) MKL) sm[K(L - z)), (A.ll) 

n=0 

where r and z are the radial distance and vertical height of the source, K (x) is the modified 
Bessel function of the second type, k v = V C /2D, and k n is the solution of the equation 
2k n cos(k n L) = — (2hT tot /D + 2k v ) sm(k n L), and c n = 1 — sm ( knL ) c ° s ( k ™ L ) _ p or ari y source 
function q(r, z, 6; E), the local observed flux is 



$|(£) = _ x2 dz rdrGUr,z,E) d6q(r, z, 9; E). (A.12) 
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